#include "pr/Pca.h"
#include "pr/PatternOperation.h"

namespace pr {

Pca::Pca(TrainingSet& trainingSet)
  : trainingSet_(trainingSet)
{
}

void Pca::computeFeatureSpace()
{
  std::size_t numFeatures;
  Pattern eigenvalues;
  std::vector<Pattern> eigenvectors;
  klt(numFeatures, eigenvalues, eigenvectors);
  featureSpace_ = FeatureSpace(numFeatures, eigenvalues, eigenvectors);
}

void Pca::transformTrainingSet()
{
//  project(mseigv);
//  trainingSet_.transform(PatternOperation::Normalize2());
}

void Pca::klt(std::size_t& numFeatures, Pattern& mseigval, std::vector<Pattern>& mseigv) const
{
//  class CenterAroundMean : public Pattern::Operation
//  {
//  public:
//    CenterAroundMean(const Pattern& mean)
//      : mean_(mean)
//    {}
//    virtual void operator()(Pattern& pattern) {
//      pattern -= mean_;
//    }
//  private:
//    Pattern mean_;
//  };
//  trainingSet_.transform(CenterAroundMean());
//
//  std::size_t patternCount = trainingSet_.patternCount();
//  std::size_t num = std::min(patternCount, trainingSet_.featureCount());
//
//  smatrix covariance(num);
//
//  if (patternCount < trainingSet_.featureCount())      // simple method (U'*U)
//  {
//    for (std::size_t i = 0; i < patternCount; ++i)
//      for (std::size_t j = i; j < patternCount; ++j)
//        covariance.utvalue(i,j) = tsetvector[i].scalarproduct(tsetvector[j]);
//  }
//  else                        // standard method (U*U')
//    covariancefromdiff(covariance);
//
//  squarematrix eigvectors(num);
//  Pattern eigvalues(num);
//
//  covariance.eig(eigvectors,eigvalues); // covariance sarà modificata
//
//  // ordino gli indici in ordine di autovalori decrescenti
//  std::vector<std::size_t> index(num);
//  for (std::size_t i = 0; i < num; ++i)
//    index[i] = i;
//  for (std::size_t i = 0; i < (num-1); ++i)
//    for (std::size_t j = i+1; j < num; ++j)
//      if (fabs(eigvalues[index[i]]) < fabs(eigvalues[index[j]]))
//      {
//        std::size_t temp = index[i];
//        index[i] = index[j];
//        index[j] = temp;
//      }
//
//  // calcolo il numero di autovettori da usare
//  util::Real sum(0.0), ro(0.0);
//  for (std::size_t i = 0; i < num; ++i)
//    sum += eigvalues[i];
//
//  std::size_t i = 0;
//  for ( ; ro/sum < precision_; ++i)
//    ro += eigvalues[ index[i] ];
//  numFeatures = i;
//  mseigv.resize(numFeatures, pvector(nfeature) );
//  mseigval.resize(numFeatures);
//
//    // creo la matrice degli autovettori più significativi
//  if (wo3)
//    h = 3; // else h=0
//
//  index.resize(nf+h);
//  if (patternCount < trainingSet_.featureCount())       // simple method
//  {
//    for (i=0; i<numFeatures; i++)         /* Matrix multiplication: O(n^3) */
//    {
//      for (j=0; j<nfeature; j++)
//        for (std::size_t k = 0; k < patternCount; ++k)                 // inline function
//          mseigv[i][j] += tsetvector[k][j]*eigvectors.ivalue(k,index[i+h]);
//      mseigv[i].normalize2();
//    }
//  }
//  else                         // standard method
//  {
//    for (i=0;i<nf; i++)
//      mseigv[i] = eigvectors.column(index[i+h]);
//  }
//  for (i=0; i<nf; i++)
//    mseigval[i] = eigvalues[ index[i+h] ];
}

} // namespace pr
